####--------------------------------------------####
# Title: Semi-parametric Selection Models for Potentially Non-ignorable 
# Attrition in Panel Studies with Refreshment Samples
# Manuscript ID PA-2013-094
# Y.SI, J.REITER AND S.HILLYGUS
####--------------------------------------------####
###--Compare semi-parametric AN model with Amelia--###
###--Plots for simulation study--###
####--------------------------------------------####
###Author: YAJUAN SI
###Latest edit date: 03/30/2014
####--------------------------------------------####
library(Amelia)
library(ggplot2)
####-------load data files from Amelia and AN model------####

load(file="output/simulation-amelia.RData") #output from Amelia
an <-read.csv('output/simumlation-an.data',header=FALSE) #output from AN model

#temparary variable name for Y2 related quantities
xvar<-c('Y21','Y22','Y23','Y24','Y25','Y11-Y21','Y12-Y22','Y13-Y23','Y14-Y24',
        'Y15-Y25','Y21-Y25','Y22-Y25','Y23-Y25','Y24-Y25')


bias_out_all<-matrix(c(apply(q_mi,2,mean)-true_q, an[,2]),length(true_q),2)
#exclude Y1
bias_out<-abs(c(bias_out_all[6:length(true_q),]))
bias_out<-data.frame(cbind(c(rep("Amelia",14),rep("Semi-AN",14)), rep(xvar,2),bias_out))
names(bias_out)<-c("type","quantity","bias")
bias_out_p<-bias_out


bias_out_p[,3]<-abs(as.numeric(as.character(bias_out_p[,3])))


bias_out_p$quantity <- as.character(bias_out_p$quantity)
bias_out_p$quantity <- factor(bias_out_p$quantity, levels=unique(bias_out_p$quantity), ordered=TRUE)

postscript("output/selection-sm-pr-bias.eps")
ggplot(data = bias_out_p, aes(x = quantity, y = bias,
                              color=type,shape=type,size=type))+
  geom_point() +
  scale_colour_manual(values=c("red","blue")) +   
  scale_size_manual(values=c(3,3))+ 
  scale_shape(solid = TRUE)+
  scale_fill_discrete("")+
  theme(
    legend.position="bottom",
    legend.direction="horizontal",
    legend.text = element_text(size = 16, face = 'bold'),
    legend.title = element_blank()
  )+
  scale_y_continuous("DIF") +
  scale_x_discrete("",                  
labels=c(expression(paste(Y[21])),expression(paste(Y[22])), expression(paste(Y[23])),expression(paste(Y[24])),
         expression(paste(Y[25])),expression(paste(Y[11-21])),expression(paste(Y[12-22])),
         expression(paste(Y[13-23])),expression(paste(Y[14-24])),
         expression(paste(Y[15-25])),expression(paste(Y[21-25])),
         expression(paste(Y[22-25])),expression(paste(Y[23-25])),
         expression(paste(Y[24-25])))) +   
  theme(
  plot.background = element_blank(),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text.x = element_text(colour="grey20",size=19,angle=0,hjust=.5,vjust=.5,face="bold"),
        axis.text.y = element_text(colour="grey20",size=15,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=0,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()


rmse_out_all<-matrix(c(rmse, an[,3]),length(true_q),2)
rmse_out<-abs(c(rmse_out_all[6:length(true_q),]))
rmse_out<-data.frame(cbind(c(rep("Amelia",14),rep("Semi-AN",14)), rep(xvar,2),rmse_out))
names(rmse_out)<-c("type","quantity","rmse")
rmse_out_p<-rmse_out

rmse_out_p[,3]<-abs(as.numeric(as.character(rmse_out_p[,3])))

rmse_out_p$quantity <- as.character(rmse_out_p$quantity)
rmse_out_p$quantity <- factor(rmse_out_p$quantity, levels=unique(rmse_out_p$quantity), ordered=TRUE)

postscript("output/selection-sm-pr-rmse.eps")
ggplot(data = rmse_out_p, aes(x = quantity, y = rmse,
                              color=type,shape=type,size=type))+
  geom_point() +
  scale_colour_manual(values=c("red","blue")) +   
  scale_size_manual(values=c(3,3))+ 
  scale_shape(solid = TRUE)+
  scale_fill_discrete("")+
  theme(
    legend.position="bottom",
    legend.direction="horizontal",
    legend.text = element_text(size = 16, face = 'bold'),
    legend.title = element_blank()
  )+
  scale_y_continuous("RMSE") +
  scale_x_discrete("",labels=c(expression(paste(Y[21])),expression(paste(Y[22])), expression(paste(Y[23])),expression(paste(Y[24])),
                               expression(paste(Y[25])),expression(paste(Y[11-21])),expression(paste(Y[12-22])),
                               expression(paste(Y[13-23])),expression(paste(Y[14-24])),
                               expression(paste(Y[15-25])),expression(paste(Y[21-25])),
                               expression(paste(Y[22-25])),expression(paste(Y[23-25])),
                               expression(paste(Y[24-25])))) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text.x = element_text(colour="grey20",size=19,angle=0,hjust=.5,vjust=.5,face="bold"),
        axis.text.y = element_text(colour="grey20",size=15,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=0,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()


cover_out_all<-matrix(c(apply(cover,2,mean), an[,1]),length(true_q),2)
cover_out<-abs(c(cover_out_all[6:length(true_q),]))
cover_out<-data.frame(cbind(c(rep("Amelia",14),rep("Semi-AN",14)), rep(xvar,2),cover_out))
names(cover_out)<-c("type","quantity","cover")
cover_out_p<-cover_out

cover_out_p[,3]<-abs(as.numeric(as.character(cover_out_p[,3])))

cover_out_p$quantity <- as.character(cover_out_p$quantity)
cover_out_p$quantity <- factor(cover_out_p$quantity, levels=unique(cover_out_p$quantity), ordered=TRUE)

postscript("output/selection-sm-pr-cr.eps")
ggplot()+
  geom_boxplot(data = cover_out_p, aes(x = type, y = cover)) +
  scale_fill_discrete("")+
  theme(
    legend.position="bottom",
    legend.direction="horizontal",
    legend.text = element_text(size = 12, face = 'bold'),
    legend.title = element_blank()
  )+
  scale_y_continuous("Coverage rates") +
  scale_x_discrete("") +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text.x = element_text(colour="grey20",size=25,angle=0,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=15,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=20,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()

#variance
t_out_all<-matrix(c(apply(t_mi,2,mean), an[,4]),length(true_q),2)

t_out<-abs(c(t_out_all[6:length(true_q),]))
t_out<-data.frame(cbind(c(rep("Amelia",14),rep("Semi-AN",14)), rep(xvar,2),t_out))
names(t_out)<-c("type","quantity","var")
t_out_p<-t_out

t_out_p[,3]<-abs(as.numeric(as.character(t_out_p[,3])))

t_out_p$quantity <- as.character(t_out_p$quantity)
t_out_p$quantity <- factor(t_out_p$quantity, levels=unique(t_out_p$quantity), ordered=TRUE)

ostscript("output/selection-sm-pr-t.eps")
ggplot(data = t_out_p, aes(x = quantity, y = var,
                              color=type,shape=type,size=type))+
  geom_point() +
  scale_colour_manual(values=c("red","blue")) +   
  scale_size_manual(values=c(3,3))+ 
  scale_shape(solid = TRUE)+
  scale_fill_discrete("")+
  theme(
    legend.position="bottom",
    legend.direction="horizontal",
    legend.text = element_text(size = 16, face = 'bold'),
    legend.title = element_blank()
  )+
  scale_y_continuous("Variance") +
  scale_x_discrete("",labels=c(expression(paste(Y[21])),expression(paste(Y[22])), expression(paste(Y[23])),expression(paste(Y[24])),
                               expression(paste(Y[25])),expression(paste(Y[11-21])),expression(paste(Y[12-22])),
                               expression(paste(Y[13-23])),expression(paste(Y[14-24])),
                               expression(paste(Y[15-25])),expression(paste(Y[21-25])),
                               expression(paste(Y[22-25])),expression(paste(Y[23-25])),
                               expression(paste(Y[24-25])))) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text.x = element_text(colour="grey20",size=19,angle=0,hjust=.5,vjust=.5,face="bold"),
        axis.text.y = element_text(colour="grey20",size=15,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=0,angle=0,hjust=.5,vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"))
dev.off()


